function [summary] = spike_sum(DFILE, voltage, isilat, abs_lat, rmp, filename, frec, dead_win, min_isi, t_left, t_right, plot_flag)
% compute information about spikes detected in the isilat table
%
% 11/28/98 P. Manis
%
% 11/12/99. P. Manis
% added spike threshold detection algorithm, following Erisir et al.
% (J. Neurophysiol. 82: 2476, 1999).
% Threshold is = V at which ([d2v/dt2] > k * sd (baseline d2v/db2)).
% where k is spike_sd_thr (as defined below).  Erisir used k=3.
% Other measurments could be made relative to this - should they?
%
% changed so that DFILE is passed.
% only DFILE argument that is neccesary is DFILE.rate,
% so things can be easily faked if necessary.

if(~exist('plot_flag', 'var'))
   plot_flag=0;
end;
if(length(min_isi) > 1)
    min_isipost = min_isi(2);
    min_isi = min_isi(1);
else
    min_isipost = 10;
end;
v = size(isilat);  % find out how many spikes are present
summary=[];
summary.spikeno = NaN; %  save position of spike in the train
summary.measthresh= NaN; % (k) = meas_pos; % is not an array element - just one for the whole thing.
summary.deadwin= NaN; % (k) = dead_win; % these don't change, but we will use just the first element later...
summary.minisi= NaN; % (k) = min_isi;
summary.minisipost = NaN;
summary.source= NaN; % (k) = rec;
summary.half_ht= NaN; % (k) = half_ht;
summary.lat= NaN; % (k) = isilat(i).lat(j); % compute spike latency relative to start of step
summary.vthr= NaN; % (k) = voltage(rec, min(idv2max)); % first point is threshold voltage.
summary.vmax= NaN; % (k) = vmax;
summary.vmin= NaN; % (k) = vmin;
summary.apht= NaN; % (k) = vmax - rmp;
summary.dvmax= NaN; % (k) = dvmax;
summary.dvmin= NaN; % (k) = -dvmin; % note sign reversal...
summary.dvvmax= NaN; % (k) = voltage(rec, idvmax); % offset is for diff function shift
summary.dvvmin= NaN; % (k) = voltage(rec, idvmin);
summary.hwa= NaN; % (k) = find_pt(voltage(rec, pbmin+pb-1:imax), half_ht, rate, 1);
summary.hwb= NaN; % (k) = find_pt(voltage(rec, imax:pe), half_ht, rate, -1);
summary.ahpdepth= NaN; % (k) = ahpdepth;
summary.ahpmin= NaN; % (k)=ahpmin;
summary.ahpmax= NaN; % (k)=ahpmax;


% --------------------------------------------------------
s_winl = t_left;
s_winr = t_right;
spike_sd_thr=5; % number of sd's above noise for dv2
if(isempty(dead_win) || isnan(dead_win))
   dead_win = 0.25;
end
RATES = DFILE.rate .* DFILE.nr_channel / 1000;
dead_time = dead_win; % msec
p_dead = floor(dead_time./RATES);
half_ht = -20; % a default half-height...

p_winl = floor(s_winl./RATES); % get "pointer" data
p_winr = floor(s_winr./RATES);
p_abs_lat = floor(abs_lat./RATES);

w = size(voltage);
v = size(isilat);
summary = [];
time=make_time(DFILE);
meas_pos = 0.5;
k = 0;

%---------------------
% do the plot... if we are set...
%
if(plot_flag >= 0)
   
   h = findobj('Tag', 'APWave');
   if(isempty(h))
      h = figure('Tag', 'APWave', 'Name', 'Action potential Waveforms', 'NumberTitle', 'off');
   end
   figure(h);
   clf; % force clean slate
   set(h, 'Backingstore', 'off');
   subplot('Position', [0.1 0.9 0.8 0.1]); cla;
   axis([0,1,0,1])
   axis('off')
   text(0.33,0.7,sprintf('%-12s R%d',filename, frec), 'Fontsize', 8);	
   hsp1=subplot('Position', [0.1 0.5 0.8 0.4]); cla; %  axis([-s_winl s_winr -80 20]);
   hsp2=subplot('Position', [0.1 0.1 0.35 0.35]); cla; % axis([-s_winl s_winr -200 200]);
   hsp3=subplot('Position', [0.6 0.1 0.35 0.35]); cla; % axis([-s_winl s_winr -200 200]);
end;
tap=[]; % time base for aps
vap=[]; % voltage for aps
dvap=[]; % derivative for aps
pbi=[]; pbe=[]; pei=[];
vthr=[]; tt0=[];
tdvmx=[]; tdvmn=[];
tmx1=[]; tmx2=[]; tmx0=[]; tahp=[];
tdvmx2=[]; tdvmn2=[];
vthr=[];
v_thr=[];

for i=1: v(2) % for each record in the latency list
   nspikes = isilat(i).nspikes; % number of isi's (should be one less than the number of spikes)
   for j = 1:nspikes % for each spike in the list
      if(j > 1 && isilat(i).isi(j-1) < min_isi) % spikes must be preceded and followed by min isi
         break;
      end
      if(j < nspikes && isilat(i).isi(j) < min_isi)
         break
      end
      if(j > 1 && j+1 < nspikes && (isilat(i).isi(j+1) < min_isipost))
          break;
      end;
      if((j+1 < nspikes) && (isilat(i).isi(j+1) < min_isipost))
         break
      end

      rec = isilat(i).source; % get trace number for spike
      if(nspikes > 1)
         t0 = isilat(i).lat(j) + abs_lat;
         pos = floor(t0/RATES(rec)); % calculate position of spike in record, don't forget the offset from the step
         s = (t0/RATES(rec)) - 1; % subtract one to correct indexing
      else
         t0 = isilat(i).lat + abs_lat;
         pos = floor(t0/RATES(rec));
         s = (t0/RATES(rec)) - 1; % subtract one to correct indexing
      end
      if(s < 1)
         s = 1;
      end;
      
      pb = floor(pos - p_winl(rec)); % attempt to use available window
      if(pb < p_abs_lat(rec))
         pb = p_abs_lat(rec)+p_dead(rec); % but never less than the absolute starting point from the driving step - and add dead time to get over the transient
      else
         pb=pb+p_dead(rec);
      end;
      pe = floor(pos + p_winr(rec));
      time2 = pb:pe;
      time2 = time2*RATES(rec) - t0;
      % compute the time where the latency was detected
      my = voltage(rec, ceil(s)) -  voltage(rec, floor(s)); % take the two adjacent points
      by = voltage(rec, floor(s)) - my * floor(s);
      v_thr = my * s + by; % this just recomputed position...
      %        disp('check')
      if(pb > 0 && pe <= w(2) && pe > pb+1) % make sure these are in the data window
         dv = gradient(voltage(rec, :),RATES(rec)); % compute dv/dt - USE whole trace
         dv2 = gradient(dv,RATES(rec)); % second derivative...
         [vmax, imax] = max(voltage(rec, pb:pe));
         if(vmax < v_thr) % bad spike detection
            break;
         end
         imaxm=imax;
         imax = imax + pb - 1;
         [vmin, imin] = min(voltage(rec, imax:pe)); % find min only on falling side
         imin = imin + imax - 1;
         % find begining of rising phase of data - min of V - if its not at first point,
         % use that new min to start the dvmax search (eliminates some artifacts)
         [nvmin, pbmin] = min(voltage(rec, pb:imax));
         beg=pbmin+pb-1; % adjust for position in voltage array
         dv2thr = spike_sd_thr*std(dv2(round((2*pe+imax)/3):pe)); % 1:floor((imax-pb+1)/3))); % over preceeding period.
         [dvmax, idvmax] = max(dv(beg:imax)); % note indices are relative to whole trace
         [dv2max, idv2max] = max(dv2(beg:imax+1)); % find(dv2(pbmin+pb-1:imax) >= dv2thr);
         dv2thr=dv2max/3; % use 1/3 threshold criterion
         [t_thr, p_thr]=find_val(dv2(beg:imax)', dv2thr, 1);
         if(~isnan(p_thr) && (p_thr < idv2max))
            p_thr=p_thr(1);
            t_thr=(p_thr+beg-1)*RATES(rec);
         else
            t_thr=NaN;
            p_thr=NaN;
         end;
         
         idvmax = idvmax + beg -1;
         [dvmin, idvmin] = min(dv(imax:pe));
         idvmin = idvmin + imax - 1;
         half_ht = rmp + meas_pos*(vmax - rmp);
         % compute ahp "depth" - compare ahp with mean of V 15-20 msec later
         [ahpmin, iahpmin] = min(voltage(rec, imax:pe));
         iahpmin = iahpmin + imax - 1;
         t_ahp = iahpmin*RATES(rec);
         te = t_ahp+15; % end of ahp post-peak measurement window
         tes = t_ahp;
         p_te = floor(te/RATES(rec));
         if(p_te > length(voltage(rec,:)))
            p_te = length(voltage(rec,:));							
         end
         p_tes = floor(tes/RATES(rec));
         [ahpmax, iahpmax] = max(voltage(rec, p_tes:p_te));
         t_ahpmax = (iahpmax+tes-1)*RATES(rec);
         ahpdepth = ahpmin - ahpmax;
         % store data to summary structure
         k=k+1;
         summary.spikeno(k)=j; % save position of spike in the train
         summary.measthresh(k) = meas_pos; % is not an array element - just one for the whole thing.
         summary.deadwin(k) = dead_win; % these don't change, but we will use just the first element later...
         summary.minisi(k) = min_isi;
         summary.minisipost(k) = min_isipost;
         summary.source(k) = rec;
         summary.lat(k) = isilat(i).lat(j); % compute spike latency relative to start of step
         if(~isnan(t_thr) && ~isnan(p_thr))
            summary.tthr(k) = (p_thr+beg-1)*RATES(rec);
            summary.vthr(k) = interp1(time(rec, :), voltage(rec, :), summary.tthr(k)); % first point is threshold voltage.
         else
            summary.tthr(k)=NaN;
            summary.vthr(k)=NaN;
         end
         summary.vmax(k) = vmax;
         summary.vmin(k) = vmin;
         summary.apht(k) = vmax - rmp;
         summary.dvmax(k) = dvmax;
         summary.dvmin(k) = -dvmin; % note sign reversal...
         summary.dvvmax(k) = voltage(rec, idvmax); % offset is for diff function shift
         summary.dvvmin(k) = voltage(rec, idvmin);
         half_ht = summary.vthr(k)+(summary.vmax(k)-summary.vthr(k))/2;
         summary.half_ht(k) = half_ht;
         summary.hwa(k) = find_pt(voltage(rec, pbmin+pb-1:imax), half_ht, RATES(rec), 1);
         summary.hwb(k) = find_pt(voltage(rec, imax:pe), half_ht, RATES(rec), -1);
         summary.ahpdepth(k) = ahpdepth;
         summary.ahpmin(k)=ahpmin;
         summary.ahpmax(k)=ahpmax;
         tdvmax = idvmax*RATES(rec);
         tdvmin = idvmin*RATES(rec);
         %				tthr = idv2max*RATES(rec);
         pbi(k)=pb;
         pei(k)=pe;
         vthr(k)=v_thr;
         tt0(k)=t0;
         tdvmx(k)=tdvmax-t0;
         tdvmn(k)=tdvmin-t0;
         tmax=imax*RATES(rec);
         tmx0(k)=tmax-t0;
         tmx1(k)=tmax-summary.hwa(k)-t0;
         tmx2(k)=tmax+summary.hwb(k)-t0;
         tahp(k)=t_ahp-t0;
         % the next section summarizes the ap waveforms on a graph, with markers for each of the measurements
         %           show_ap(summary, k, voltage);
         tdvmx2(k)=tdvmax-t0;
         tdvmn2(k)=tdvmin-t0;
         tthr2(k)=summary.tthr(k)-t0+RATES(rec);
         vdv2thr(k)=dv2thr;
         
      end % of if statement - spikes in the window
   end % of loop on each spike in the
end

if(k == 0)
   return;
end;
% quick loop over k to get the data back and fill the arrays.
pmax=max(pei-pbi)+1; % this is the max size of the array we need.
for i=1:k
   rec = summary.source(i);
   time2=pbi(i):pei(i);
   time2=(time2*RATES(rec)) - tt0(i);
   n=size(time2);
   tap(:,i) = pad(time2, pmax)';
   xmn = min(tap(:,i));
   xmx = max(tap(:,i));
   vtemp  = voltage(rec, pbi(i):pei(i));
   ve(i)=voltage(rec,pei(i));
   vap(:,i)=pad(vtemp,pmax)';
   vdvg  = gradient(vap(:,i),RATES(rec));
   vdv(:,i) = pad(vdvg, pmax);
   
   vdv2g=gradient(vdv(:,i),RATES(rec));
   vdv2(:,i)=pad(vdv2g, pmax);
   %         disp(sprintf('i: %3d r: %3d t0: %6.1f t1: %6.1f', i, rec, xmn, xmx))
end
xmin = min(min(tap));
xmax = max(max(tap));
xl = [xmin xmax];
% now we can graph the data FAST...
if(plot_flag >= 0)
   subplot(hsp1);
   axis auto
   if(~isempty(vap))
      plot(tap, vap, 'Parent', hsp1, 'color', 'black', 'LineStyle', '-');
      xlim(xl);
      text(-1, v_thr, 'lat\rightarrow',  'HorizontalAlignment', 'right', 'color', 'red'); % latency detection position
      text(tthr2, summary.vthr, '\rightarrow',  'HorizontalAlignment', 'right', 'color', 'blue'); % threshold detection position
      text(tmx0, summary.vmax, '\downarrow',  'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'color', 'blue');
      text(tmx1-1, summary.half_ht, '\leftrightarrow', 'Parent', hsp1, ...
         'HorizontalAlignment', 'right', 'color', 'magenta');
      text(tmx2, summary.half_ht,  '\uparrow', ...
         'Parent', hsp1 ,'HorizontalAlignment', 'center', 'VerticalAlignment', 'top', 'color', 'magenta');
      %         else
      %           text(max(time2), ve, 'no HWB', 'color', 'r', 'VerticalAlignment', 'middle', 'HorizontalAlignment', 'left', 'Fontsize', 7);
      %       end
      text(tahp, summary.ahpmin, '\uparrow', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'top', 'color', 'g');
      text(tahp, summary.ahpmax, '\downarrow',  'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'color', 'g');
      text(tdvmx, summary.dvvmax, 'X',  'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'color', 'r');
      text(tdvmn, summary.dvvmin,  'X', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'color', 'r');
   else
      text(0.5, 0.5, 'No Spikes', 'HorizontalAlignment', 'center');
   end;
   
   
   subplot(hsp2);
   axis auto
   if(~isempty(vdv))
      plot(tap, vdv, 'Parent', hsp2, ...
         'color', 'black', 'LineStyle', '-');
      xlim([xl(1) xl(2)/4]);
      text(tdvmx, summary.dvmax,  '\downarrow', 'Parent', ...
         hsp2, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'color', 'red');
      text(tdvmn, -summary.dvmin, '\uparrow',  'Parent', ...
         hsp2, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'top', 'color', 'red');
   end;
   subplot(hsp3);
   axis auto
   if(~isempty(vdv2))
      plot(tap, vdv2, 'Parent', hsp3, ...
         'color', 'black', 'LineStyle', '-');
      xlim([xl(1) xl(2)/4]);
      text(tthr2, vdv2thr, '\rightarrow',  'HorizontalAlignment', 'right', 'color', 'green'); % threshold detection position
   end;
   
   drawnow;
   
   
   %print -dljet3;
end;

return;

function [newt, p] = find_pt(v, t, rate, dir)
% This function finds the first point to the right in v that is >= t if dir = 1
% or <= t if dir = -1, and interpolates to find the time at which t would occur
% the t value is returned in newt...

newt=NaN;
p = NaN;	 % optional return argument for found position (not integer!)

npts = length(v);
if(npts < 2)
   return;
end

if(dir>0)
   v0 = max(v, t); % identify points > t
   v1 = find(v0 > t); % make it a short list
   if(isempty(v1))
      return;
   end
   v2=v1(1);   % get just the first point
   if(v2 > 1)
      m = v(v2) - v(v2-1); % compute the slope
   else
      m = v(2)-v(1);
   end
   b = v(v2) - m*v2; % compute the intercept (note, using points at time units...)
   
   pt = ((t - b) / m); % correct for indexing - but is calculated relative to start of the data passed (e.g., pb)
   newt = (npts-pt)*rate;
else
   v0 = max(v, t); % identify points in v that are > t
   v1 = find(v0 == t); % make it a short list - of the points that were <= t!
   if(isempty(v1)) % bad points...
      return;
   end
   v2=v1(1);   % get the first point that fell to or below the threshold
   if(v2 < length(v))
      m = v(v2+1) - v(v2); % compute the slope
   else
      m = v(v2)-v(v2-1); % one way or the other...
   end
   b = v(v2) - m*v2; % compute the intercept (note, using points at time units...)
   newt = (((t - b) / m)*rate)-rate; % correct for indexing
end
if(nargout > 1)
   p = (t - b)/m;
end;
return;


function [v, p] = find_val(A, b, slope)
% find point nearest b in A, in region with appropriate slope
% slope value is any positive negative or 0 value

v=[];
p=[];
% first mask A off so that we look only at regions with the proper slope value
dv=gradient(A);
if(slope > 0)
   a=find(dv<0); % mask off points with negative slope
   mask = -Inf;
else
   a=find(dv>0); % maks off points with positive slope
   mask = Inf;
end;

A(a)=mask; % masks those points off.
if(any(A~=mask))
   [v, p] = findnear(A, b);
end;
return;
